%     VIRTUAL SORIC - программа для получения эффективных свойств заполнителей SORIC
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [NODES, ELEMENTS, RVE_Coordinates, Node_Merge_Tolerance] = Generate_Mesh(Points, N_Elements_Edge, RVE_Thickness, N_Elements_High, Shell_Layer, Shell_Layer_Thickness, RVE_Size, Element_Order, RVE_Array_1, RVE_Array_2, RVE_Array_3)

global Merge_Tolerance_Multiplier;
global Element_Critical_Ratio;

global ErrCode;

ELEMENTS = [];
RVE_Coordinates = [];
Node_Merge_Tolerance = [];

fprintf (1, 'Построение сетки...                                     ');

% Построение узлов на плоскости
Nodes_2D = Construct_Nodes_2D(Points, N_Elements_Edge);

% Построение узлов в объеме
NODES = Construct_Nodes_3D(Nodes_2D, N_Elements_Edge, RVE_Thickness, N_Elements_High, Shell_Layer, Shell_Layer_Thickness);

% Построение элементов C3D8
EL_C3D8 = Make_EL_C3D8(NODES, N_Elements_Edge, N_Elements_High, Shell_Layer);

% Построение элементов C3D6
EL_C3D6 = Make_EL_C3D6(NODES, N_Elements_Edge, N_Elements_High, Shell_Layer);

% Минимальная длина ребра элемента
Min_Element_Edge = min([norm(Points(27,:) - Points(15,:)) / N_Elements_Edge, ...
                        norm(Points(15,:) - Points(16,:)) / N_Elements_Edge, ...
                        norm(Points(16,:) - Points(28,:)) / N_Elements_Edge, ...
                        norm(Points(28,:) - Points(27,:)) / N_Elements_Edge, ...
                        norm(Points(16,:) - Points(6,:)), ...
                        norm(Points(6,:) - Points(17,:)), ...
                        norm(Points(17,:) - Points(16,:)), ...
                        norm(Points(15,:) - Points(14,:)), ...
                        norm(Points(14,:) - Points(7,:)), ...
                        norm(Points(7,:) - Points(15,:))]);
                    
if (strcmp(Shell_Layer, 'TRUE'))
    Min_Vertical_Edge = min(Shell_Layer_Thickness, (RVE_Thickness - 2*Shell_Layer_Thickness) / N_Elements_High);
else
    Min_Vertical_Edge = RVE_Thickness / N_Elements_High;
end

if (Min_Vertical_Edge < Min_Element_Edge)
    Min_Element_Edge = Min_Vertical_Edge;
end

% Контроль соотношения сторон элементов
Element_Ratio = (norm(Points(27,:) - Points(15,:)) / N_Elements_Edge) / ...
                (norm(Points(27,:) - Points(28,:)) / N_Elements_Edge);
            
if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end
            
Element_Ratio = (norm(Points(16,:) - Points(15,:)) / N_Elements_Edge) / ...
                (norm(Points(28,:) - Points(16,:)) / N_Elements_Edge);

if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end

Element_Ratio = (norm(Points(28,:) - Points(16,:)) / N_Elements_Edge) / ...
                norm(Points(28,:) - Points(29,:));
            
if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end
            
Element_Ratio = norm(Points(28,:) - Points(29,:)) / ...
                (norm(Points(28,:) - Points(16,:)) / N_Elements_Edge);

if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end

Element_Ratio = (norm(Points(16,:) - Points(15,:)) / N_Elements_Edge) / ...
                norm(Points(16,:) - Points(6,:));
            
if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end
            
Element_Ratio = norm(Points(16,:) - Points(6,:)) / ...
                (norm(Points(16,:) - Points(15,:)) / N_Elements_Edge);

if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end

Element_Ratio = (norm(Points(14,:) - Points(26,:)) / N_Elements_Edge) / ...
                norm(Points(14,:) - Points(7,:));
            
if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end
            
Element_Ratio = norm(Points(14,:) - Points(7,:)) / ...
                (norm(Points(14,:) - Points(26,:)) / N_Elements_Edge);

if (Element_Ratio > Element_Critical_Ratio)
    fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое соотношение сторон конечных элементов. Измените параметры построения сетки.\n');
    ErrCode = -1;
    return;
end

% Допуск для слиянии узлов
Node_Merge_Tolerance = Min_Element_Edge * Merge_Tolerance_Multiplier;

% Создание единого массива элементов
ELEMENTS = zeros(size(EL_C3D8,1) + size(EL_C3D6,1), 12);

ELEMENTS(1:size(EL_C3D8,1), 1) = EL_C3D8(:,1);
ELEMENTS(size(EL_C3D8,1)+1:size(ELEMENTS,1), 1) = EL_C3D6(:,1);

ELEMENTS(1:size(ELEMENTS,1), 2) = 3;

ELEMENTS(1:size(EL_C3D8,1), 3) = 8;
ELEMENTS(size(EL_C3D8,1)+1:size(ELEMENTS,1), 3) = 6;

ELEMENTS(1:size(EL_C3D8,1), 4) = EL_C3D8(:,2);
ELEMENTS(size(EL_C3D8,1)+1:size(ELEMENTS,1), 4) = EL_C3D6(:,2);

ELEMENTS(1:size(EL_C3D8,1), 5:12) = EL_C3D8(:,3:10);
ELEMENTS(size(EL_C3D8,1)+1:size(ELEMENTS,1), 5:10) = EL_C3D6(:,3:8);

% Габариты ячейки
RVE_Coordinates = zeros(3,2);

RVE_Coordinates(1,1) = 0;
RVE_Coordinates(1,2) = RVE_Size(1);
RVE_Coordinates(2,1) = 0;
RVE_Coordinates(2,2) = RVE_Size(2);
RVE_Coordinates(3,1) = 0;
RVE_Coordinates(3,2) = RVE_Size(3);

fprintf (1, 'OK!\n');

% Конвертация элементов в элементы второго порядка
if (strcmp(Element_Order, 'QUADRATIC'))
    
    fprintf (1, 'Конвертация элементов в квадратичные...                 ');
    
    % Уменьшение допуска для слияния узлов
    Node_Merge_Tolerance = Node_Merge_Tolerance/2;
    
    [NODES, ELEMENTS] = Convert_Mesh_To_Quad(NODES, ELEMENTS, Node_Merge_Tolerance);
    
    fprintf (1, 'OK!\n');
    
end

% Создание массива ячеек
if ((RVE_Array_1 > 1) || (RVE_Array_2 > 1) || (RVE_Array_3 > 1))
    
    fprintf (1, 'Создание массива ячеек...                               ');
    
    [NODES, ELEMENTS, RVE_Coordinates] = Create_RVE_Array(NODES, ELEMENTS, Node_Merge_Tolerance, RVE_Array_1, RVE_Array_2, RVE_Array_3, RVE_Size, RVE_Coordinates);
    
    fprintf (1, 'OK!\n');
    
end
